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Abstract 

Point-based  surface  processing  has  developed  into  an  at¬ 
tractive  alternative  to  mesh-based  processing  techniques 
for  a  number  of  geometric  modeling  applications.  By 
working  with  point  cloud  data  directly,  any  processing 
is  based  on  the  given  raw  data  and  its  underlying  geome¬ 
try  rather  than  any  arbitrary  intermediate  representations 
and  generally  artificial  connectivity  relations.  In  this  pa¬ 
per,  we  introduce  the  notion  of  meshless,  or  point  cloud, 
subdivision  by  extending  concepts  from  recursive  mesh 
subdivision  to  the  point  cloud  case.  We  are  primarily 
concerned  with  showing  the  conceptual  viability  of  this 
idea  and  propose  a  first  geometric  meshless  subdivision 
framework.  By  replacing  the  role  of  mesh  connectiv¬ 
ity  by  intrinsic  point  proximity  information  and  devis¬ 
ing  a  meshless  geodesic  subdivision  operator,  we  avoid 
the  costly  surface  reconstruction,  simplification  and  po¬ 
tential  remeshing  preprocessing  steps  typically  required 
for  supporting  mesh-based  subdivision,  steps  which  are 
in  general  not  directly  related  to  the  underlying  object 
geometry.  Furthermore,  the  maintenance  of  any  global 
combinatorial  data  structure  such  as  a  mesh  connectivity 
graph  is  not  required.  This  property  also  makes  our  ap¬ 
proach  relatively  easily  extensible  to  the  processing  of 
point-based  representations  of  higher-dimensional  sur- 

*  The  author  performed  this  work  whilst  visiting  the  University  of 
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faces.  Apart  from  introducing  the  idea  of  meshless 
subdivision,  our  main  contributions  are,  firstly,  a  first 
meshless  geodesic  subdivision  operator.  Secondly,  we 
present  a  new  method  for  the  computation  of  geodesic 
weighted  averages  on  manifold  surfaces,  which  are  at 
the  heart  of  our  point  cloud  subdivision  framework. 

1  Introduction 

Point-based  editing,  free-form  and  multiresolution  mod¬ 
elling  and  visualisation  have  developed  into  viable  alter¬ 
natives  to  mesh-based  processing  methods  [ABCCT03, 
PGK02,  PKKG03,  ZPvBGOl,  ZPKG02],  With  nu¬ 
merous  applications  in  medical  imaging,  reverse  engi¬ 
neering,  cultural  heritage,  geoscience,  etc.  acquiring 
point  sets  of  significant  density  from  three  or  higher¬ 
dimensional  manifold  surfaces,  it  is  particularly  attrac¬ 
tive  to  be  able  to  work  with  point-sampled  geometry  di¬ 
rectly. 

In  the  context  of  mesh-based  surface  processing,  sub¬ 
division  has  developed  into  a  powerful  and  widely- 
used  tool  for  the  free-form  design  and  representation 
of  smooth  surfaces.  These  schemes  use  a  subdivision 
operator  which  is  recursively  applied  to  a  coarse  base 
mesh.  As  a  result,  a  sequence  of  refined  meshes  is  ob¬ 
tained  which  quickly  converges  to  a  smooth  limit  sur¬ 
face.  For  many  applications,  the  initial  base  mesh  is 
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obtained  by  costly  polygonisation  of  point-sampled  ge¬ 
ometry.  Highly  dense  point  sets  thereby  translate  into 
excessively  dense  polygonal  models  with  arbitrary  con¬ 
nectivity.  This  is  generally  followed  by  another  costly 
mesh  simplification  step  and  often  the  need  for  remesh¬ 
ing  ([ACSD*03]  and  the  references  therein)  to  enforce 
subdivision  connectivity  of  the  base  mesh.  Ease  of  ma¬ 
nipulation  tends  to  be  further  affected  by  the  inherent 
overhead  of  mesh  data  structure  maintenance  and  traver¬ 
sal. 

This  paper  presents  a  first  framework  for  meshless,  or 
point  cloud,  subdivision.  That  is,  we  propose  to  avoid 
the  consideration  of  mesh  connectivity  graphs  and  in¬ 
stead  to  work  with  the  given  point-sampled  geometry 
directly  and  intrinsically  throughout.  We  undertake  a 
first  step  towards  such  an  extension  of  recursive  mesh 
subdivision  to  point  clouds  by  introducing  a  meshless 
geodesic  subdivision  operator  reminiscent  of  widely- 
used  mesh  subdivision  operators.  More  specifically,  we 
replace  mesh  connectivity  by  intrinsic  proximity  infor¬ 
mation  and  introduce  meshless  subdivision  rules  using 
geodesic  centroids  of  local  intrinsic  point  neighbour¬ 
hoods.  Although  some  of  these  geometric  operations 
may  be  approximated  in  an  Euclidean  context  when 
working  with  large  sampling  densities,  regular  meshes 
and  very  local  subdivision  rules,  working  with  the  point- 
sampled  geometry  directly  and  intrinsically  has  the  ben¬ 
efit  of  not  requiring  any  non-geometric  preprocessing 
steps.  Furthermore,  an  intrinsic  point  cloud  subdivision 
approach  is  more  generally  applicable  in  that  geodesic 
averaging  rules  can  be  non-local  and  do  not  vary  with 
the  particular  type  of  data  (mesh)  representation  used. 

The  focus  of  this  paper  is  on  highlighting  the  concep¬ 
tual  viability  of  meshless  subdivision  by  putting  forward 
a  first  suitable  framework.  Our  ideas  on  the  theoretical 
analysis  of  meshless  geodesic  subdivision  schemes  will 
be  reported  elsewhere. 

Both  our  intrinsic  point  cloud  simplification  method 
for  the  determination  of  a  base  point  set  and  our 
geodesic  centroid  technique  are  based  upon  the  intrinsic 
distance  mapping  algorithm  for  point  clouds  presented 
in  [MS03].  Following  a  brief  overview  of  related  work 
in  Section  2,  we  therefore  summarise  the  main  aspects 
of  this  technique  in  Section  3.  Section  4  presents  the 
various  elements  of  our  intrinsic  meshless  subdivision 
framework.  Section  5  gives  experimental  results.  In  the 
concluding  Section  6,  we  briefly  remark  on  our  idea  for 
the  theoretical  analysis  of  intrinsic  meshless  subdivision 
schemes. 


2  Related  work 

We  start  with  providing  a  brief  summary  of  the  ideas 
underpinning  mesh-based  subdivision  of  surfaces  in  R3. 
The  overview  is  motivational  in  nature,  for  a  thorough 


formal  treatment  of  mesh  subdivision,  see  [DL02], 

Following  the  notation  of  [DF02],  surface  subdivi¬ 
sion  schemes  consist  of  a  subdivision  operator  S  recur¬ 
sively  applied  to  control  nets  Ni  =  NiV1  ,El ,Fl)  of  ar¬ 
bitrary  topology,  with  /  £  Zq  denoting  the  subdivision 
level,  V  representing  a  set  of  control  vertices  in  R3  and 
E  and  F  describing  the  topological  relations  in  the  form 
of  edges  and  faces  respectively.  That  is,  the  iterative  ap¬ 
plication  of  a  subdivision  scheme  generates  a  sequence 

N,+i  =  SN,. 

Define  the  surface  X  C  R3  as  the  limit  surface  of  a 
scheme  given  that 

lim  dniy1  ,X)  =  0, 

l— *oo 

with 

dn(X,Y)  =  max{sup  inf  \\x  —  y||2,sup  inf  ||x  —  v||i} 

'  xexy£Y 

representing  the  Euclidean  Hausdorff  distance  between 
the  closed  subsets  X,Y  C  R3.  Starting  with  the  coarse 
base  net  No,  at  each  iteration,  new  control  vertices  are 
inserted  and  connected  according  to  the  scheme’s  re¬ 
finement  rule.  Control  vertices  are  then  re-positioned 
following  the  operator’s  geometric  averaging  rule.  Both 
the  refinement  and  the  geometric  rule  give  the  position 
of  control  vertices  in  N/+  \  in  the  form  of  weighted  aver¬ 
ages  of  topologically  neighbouring  vertices  in  A/.  The 
careful  choice  of  these  rules  in  relation  to  the  valency 
of  the  control  vertex  under  consideration  guarantees  the 
convergence  of  the  scheme  to  a  limit  surface  of  prov¬ 
able  continuity.  Note  that  for  the  subdivision  schemes 
in  the  literature  the  convergence  is  in  each  component 
and  in  the  uniform  norm.  Convergence  in  the  Hausdorff 
distance  is  not  yet  achieved. 

Not  every  existing  mesh  subdivision  operator  allows 
for  such  a  simple  distinction  between  a  topological  re¬ 
finement  and  a  geometric  averaging  rule  applicable  at 
all  points.  However,  those  that  do  allow  for  this  kind  of 
distinction,  include  the  most  widely-used  schemes.  For 
example,  Foop  [Foo87]  subdivision  for  triangular  con¬ 
trol  nets  may  be  cast  in  this  form.  In  the  case  of  Foop 
subdivision,  the  refinement  rule  consists  of  the  insertion 
of  a  new  point  at  the  midpoint  of  every  edge,  while  the 
geometric  rule,  applicable  at  all  points  in  the  new  mesh, 
smoothes  the  locations  of  the  points  by  weighted  averag¬ 
ing  of  their  topological  neighbours  in  the  refined  mesh. 

In  this  paper,  we  propose  to  replace  this  use  of  mesh 
connectivity  by  intrinsic  proximity  information  and  for¬ 
mulate  meshless  refinement  and  geometric  averaging 
rules  in  the  form  of  weighted  geodesic  centroids  of  local 
neighbourhoods.  As  a  result,  we  obtain  a  meshless  sub¬ 
division  scheme  reminiscent  of  mesh-based  subdivision 
operators. 
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Figure  1:  Intrinsic  distance  mapping  using  Fast  March¬ 
ing  for  point  clouds  operates  in  an  offset  band  consisting 
of  the  union  of  balls  B(pi,r)  centred  at  (black)  points  p, 
of  the  surface  M  (left).  Only  those  (blue)  grid  points 
falling  inside  the  offset  band  are  considered  during  pro¬ 
cessing.  Cross-sectional  view  of  a  constant  radius  offset 
band  for  the  Michelangelo  Youthful  data  set  (right). 


To  the  best  of  our  knowledge,  there  exists  only 
one  previous  piece  of  work  which  touches  upon  this 
notion  of  meshless  subdivision.  In  Fleishman  et 
al.  [FCOAS03],  the  authors  generate  progressive  levels- 
of-detail  of  point  clouds  by  transferring  the  mesh-based 
idea  of  subdivision  displacement  maps  to  the  point 
cloud  case.  They  devise  a  point  cloud  simplification 
method  for  the  generation  of  a  base  point  set  and  present 
both  a  projection  and  a  local,  uniform  upsampling  oper¬ 
ator  with  the  help  of  local  surface  reconstruction  using 
Moving  Least  Squares  [ABCO*03],  By  contrast,  we  are 
not  aiming  to  mimic  the  principle  of  mesh-based  sub¬ 
division  displacement  mapping  for  point-based  surface 
representations  but  rather  are  interested  in  transferring 
the  idea  of  mesh  subdivision  to  the  point  cloud  case. 

In  the  following,  we  propose  a  first  intrinsic  frame¬ 
work  for  the  meshless  subdivision  of  a  point  cloud  P ) 
to  a  refined  point  cloud  P/+i,  with  /  defined  as  above. 
We  start  our  discussion  with  a  summary  of  the  intrinsic 
distance  mapping  algorithm  [MS03]  at  the  heart  of  this 
intrinsic  framework.  This  summary  is  followed  by  the 
presentation  of  the  intrinsic  meshless  subdivision  frame¬ 
work  itself. 


3  Intrinsic  distance  mapping 
across  point  clouds 

We  summarise  the  extension  of  the  well-known  origi¬ 
nal  Fast  Marching  level  set  method  [HPCD96,  Set99, 
Tsi95]  to  the  case  of  surfaces  in  point  cloud  form 
as  introduced  in  [MS03].  Our  review  is  necessarily 
terse,  presenting  just  the  key  results.  For  full  details, 
see  [MS03], 

Let  P  =  {pitPIt  ■  ■  iPn]  denote  a  set  of  points,  or 
point  cloud,  acquired  from  a  smooth,  compact  manifold 


surface  M  in  m  >  3  dimensions.  Define  the  /--offset  £lp 
as  the  union  of  Euclidean  balls  with  radius  r  centred  at 
points  pi  €  P 

n 

Q.'p  :=  U  B(pi,r)  =  {x  e  K'"  :  d(pi,x)  <  r}, 
i=i 

where  d(.,.)  denotes  the  Euclidean  distance  function 
(Figure  1).  To  approximate  the  weighted  intrinsic  dis¬ 
tance  map  originating  from  a  source  point  q  £  M  on  M, 
Memoli  and  Sapiro  [MS03]  suggest  computing  the  Eu¬ 
clidean  distance  map  in  £l'P.  That  is 

\VMTM{p)\=F{p),  (1) 

for  p  £  M  and  with  boundary  condition  TM(q)  =  0  is 
approximated  by 


|V7nr  (/?)|  =  F(p),  (2) 

for  p  £  ClrP  and  boundary  condition  Taj,  (q)  =  0.  F  rep¬ 
resents  the  (smooth)  extension  of  the  propagation  speed 
F  on  M  into  £l'p.  Tip)  denotes  the  arrival  time  at  p 
of  the  front  originating  from  q  and  Vm  and  V  repre¬ 
sent  the  intrinsic  and  the  Euclidean  gradient  operators 
respectively.  The  problem  of  computing  an  intrinsic  dis¬ 
tance  map  is  therefore  transformed  into  the  problem  of 
computing  an  extrinsic  (Euclidean)  distance  map  in  the 
tubular  neighbourhood  £l'p  around  the  surface,  i.e.  in  an 
Euclidean  manifold  with  boundary.  In  [MS03],  the  au¬ 
thors  prove  uniform  probabilistic  convergence  between 
these  two  distance  maps  and  thus  show  that  the  approxi¬ 
mation  error  between  the  intrinsic  and  extrinsic  distance 
maps  is  of  the  same  theoretical  order  as  that  of  the  Fast 
Marching  algorithm  for  both  noise-free  and  noisy  point 
clouds  (assumming  noise  to  be  bounded  from  above  by 
r ).  The  Fast  Marching  method  can  therefore  be  used  to 
approximate  the  solution  to  (2)  in  a  computationally  op¬ 
timal  manner  and  without  the  need  for  any  prior  surface 
reconstruction  by  only  slightly  modifying  the  original 
Cartesian  Fast  Marching  technique  to  deal  with  bounded 
spaces.  The  relatively  straightforward  implementation 
of  this  technique  consists  of,  firstly,  computing  the  offset 
band  Q.p,  followed  by  the  application  of  standard  Carte¬ 
sian  Fast  Marching  restricted  to  £l'p.  For  more  imple- 
mentational  details,  see  Memoli  and  Sapiro  [MS03]. 

This  approach  underpins  both  the  point  cloud  sim¬ 
plification  technique  for  the  generation  of  a  base  point 
set  and  our  geodesic  computation  algorithm  presented 
as  parts  of  our  meshless  subdivision  framework  in  the 
following  section. 


4  Intrinsic  meshless  subdivision 

We  start  with  a  brief  summary  of  the  intrinsic  point  set 
simplification  method  utilised  for  the  computation  of  a 
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base  point  set  Pq.  This  is  followed  by  the  considera¬ 
tion  of  a  suitable  intrinsic  neighbourhood  concept.  The 
section  concludes  with  the  presentation  of  the  meshless 
subdivision  operator  and  our  method  for  the  computa¬ 
tion  of  geodesic  centroids. 

4.1  Intrinsic  point  cloud  simplification 

Depending  on  the  method  used  for  the  acquisition  of 
a  point  cloud  representation  of  a  surface,  the  resulting 
point-sampled  geometry  might  be  extremely  dense.  We 
simplify  any  such  input  point  cloud  P  to  a  coarser  base 
point  set  Pq  in  a  preprocessing  step.  This  operation  is 
to  be  performed  subject  to  a  minimum  density  condition 
to  support  the  meaningful  computation  of  geodesic  cen¬ 
troids  across  Pq.  Note  in  this  context  that  in  contrast  to 
mesh  subdivision  preprocessing,  this  simplification  step 
is  fully  geometric  in  nature. 

As  regards  the  particular  simplification  technique 
used,  Pauly  et  al.  [PGK02]  introduce  a  number  of 
point  cloud  simplification  methods  by  adapting  various 
widely  used  mesh  simplification  techniques  to  the  point 
cloud  setting.  Although  their  particle  simulation-based 
simplification  method  represents  an  interesting  alterna¬ 
tive,  we  opt  for  the  farthest  point  simplification  scheme 
presented  in  [MD03]  instead.  This  choice  is  due  to 
the  method’s  superior  efficiency,  its  purely  intrinsic  na¬ 
ture,  its  simple  density  control  and  its  close  relationship 
with  a  useful  (intrinsic)  local  proximity  concept  (Sec¬ 
tion  4.2). 

The  point  cloud  simplification  method  [MD03]  ex¬ 
ploits  the  observation  that  the  progressive  farthest  point 
sampling  of  a  point  cloud  P  may  be  implemented 
by  incremental  intrinsic  (bounded)  Voronoi  diagram, 
BVD(P),  computation  [MD03,  OBSOO].  This  incremen¬ 
tal  computation  is  performed  following  an  expanding 
waves  view:  In  analogy  to  the  dropping  of  pebbles  in 
still  water,  circular  fronts  move  across  the  surface  from 
the  point  of  impact.  The  locations  where  wave  fronts 
meet  or  leave  the  domain  define  the  intrinsic  bounded 
Voronoi  diagram  of  the  points  of  impact.  See  Figure  2 
for  a  triangular  mesh-based  example  produced  using 
public  domain  software  [PC03],  This  wave  propaga¬ 
tion  is  discretised  and  simulated  accurately  by  solving 
the  Eikonal  equation  (1)  using  the  extended  Fast  March¬ 
ing  method  for  intrinsic  distance  mapping  across  point 
clouds  [MS03]  summarised  in  the  previous  section.  By 
letting  the  user  control  the  radius  of  the  largest  empty 
circle  in  the  domain  of  the  simplified  point  set,  this 
method  provides  guaranteed  bounds  on  the  minimum 
and  maximum  distance  between  neighbouring  sample 
points  (Figure  3).  It  thus  allows  for  the  relatively  simple 
control  of  a  guaranteed  density  of  the  simplified  point 
set  [MD03],  We  exploit  this  property  to  simplify  any 
highly  dense  input  P  to  a  base  point  set  Pq  still  suffi¬ 
ciently  dense  to  support  meaningful  geodesic  centroid 


Figure  2:  Wave  propagation  for  the  incremental  com¬ 
putation  of  discrete  intrinsic  bounded  Voronoi  diagrams 
and  thus  progressive  intrinsic  (red)  farthest  sample  point 
localisation  of  21,  23  and  25  sample  sites  on  a  triangu¬ 
lated  surface  (from  left  to  right). 

computations. 

The  use  of  this  simplification  technique  also  yields 
a  discrete  intrinsic  Voronoi  diagram  of  the  simplified 
point  set.  As  discussed  next,  the  availability  of  this  di¬ 
agram  suggests  the  collecting  of  local  proximity  infor¬ 
mation  using  a  Voronoi  diagram-based  neighbourhood 
concept. 

4.2  Intrinsic  proximity  information 

Subdivision  schemes  incorporate  refinement  and  geo¬ 
metric  averaging  rules  in  the  form  of  weighted  averages 
of  local  neighbourhoods.  Whilst  in  a  surface  mesh  con¬ 
text  local  neighbourhoods  follow  from  mesh  connectiv¬ 
ity,  in  the  meshless  case  this  connectivity  information  is 
replaced  by  local  proximity  information.  It  is  interesting 
to  note  in  this  context  how  information  which  is  gener¬ 
ally  not  related  to  the  geometry  of  the  problem  such  as 
mesh  connectivity  plays  such  an  important  role  for  mesh 
subdivision.  For  example,  the  same  object  when  repre¬ 
sented  with  different  types  of  connectivity  requires  the 
application  of  different  mesh  subdivision  schemes  al¬ 
though  it  is  geometrically  unchanged.  This  importance 
has  contributed  to  the  substantial  research  efforts  in  the 
area  of  remeshing.  The  role  of  mesh  connectivity  in 
mesh  subdivision  also  explains  the  existence  of  numer¬ 
ous  subdivision  operators  restricted  in  applicability  to  a 
particular  type  of  mesh  only.  By  contrast,  as  described 
in  the  following,  in  the  case  of  our  meshless  subdivision 
operator,  point  cloud  proximity  is  determined  intrinsi¬ 
cally  with  the  subdivision  operator  purely  formulated  in 
terms  of  local  proximity. 

To  allow  for  any  irregularity  in  the  base  point  set  Pq, 
we  favour  the  use  of  a  neighbourhood  concept  ensuring 
a  spherical  distribution  of  neighbours  g,  around  the  point 
p  under  consideration.  Simple  ball  or  k  nearest  neigh¬ 
bourhoods  fail  to  guarantee  such  a  neighbour  distribu¬ 
tion  [FR01],  Neighbourhood  concepts  meeting  this  re¬ 
quirement  include  Linsen’s  [LinOl]  enhanced  k  nearest 
neighbourhood,  which  enforces  a  maximum  angle  be¬ 
tween  successive  neighbours  around  p,  and  Floater  and 
Reimer’s  [FR01]  local  Delaunay  neighbourhood.  Since 
our  simplification  algorithm  provides  a  (discrete)  intrin- 
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Figure  3:  The  intrinsic  point  cloud  simplification  al¬ 
gorithm  [MD03]  incrementally  computes  a  geodesic 
Voronoi  diagram  across  the  point-sampled  geometry.  As 
illustrated  here  for  the  planar  case,  it  guarantees  a  user- 
controlled  density  by  letting  the  user  set  the  radius  p  of 
the  largest  empty  circle  in  the  domain  of  the  simplified 
point  set.  The  availability  of  the  geodesic  Voronoi  dia¬ 
gram  supports  the  immediate  determination  of  intrinsic 
proximity  information  in  the  form  of  natural  neighbour¬ 
hoods  (dashed  lines).  Note  the  spherical  distribution  of 
these  neighbours  all  around  the  input  point  under  con¬ 
sideration  for  the  kind  of  irregular  uniformity  resulting 
from  the  enforcement  of  p. 


sic  Voronoi  diagram  of  Pq,  we  collect  local  proximity 
information  by  considering  the  set  of  intrinsic  Voronoi, 
or  natural,  neighbours  Np  instead, 

Np  =  {q:  p  and  q  are  neighbours  in  BVD(Pi) }, 

for  p,  q  £  Pi,  p  /  q.  BVD(Pi),  and  thus  the  neighbour¬ 
hood  information,  may  be  maintained  by  updating  the 
diagram  with  any  points  inserted  into  P/  thereby  ob¬ 
taining  BVD(Pi+ 1)  (Figure  4).  Once  a  refined  point  set 
has  reached  a  certain  density,  the  natural  neighbourhood 
information  may  be  replaced  by  the  k  neighbours  in¬ 
trinsically  nearest  to  p.  Note  that  although  the  natural 
neighbourhood  concept  does  not  guarantee  a  spherical 
neighbour  distribution  when  dealing  with  extreme  irreg¬ 
ularity,  as  indicated  in  Figure  3,  it  performs  well  with 
the  uniform  irregular  density  guaranteed  by  the  intrin¬ 
sic  point  cloud  simplification  algorithm  discussed  in  the 
previous  section. 

This  neighbourhood  information  is  exploited  for  the 
computation  of  local  weighted  centroids  as  part  of  our 
meshless  subdivision  scheme  presented  next. 

4.3  An  intrinsic  subdivision  operator  for 
point  clouds 

Mesh-based  subdivision  uses  rules  based  on  weighted 
averages  of  local  neighbourhoods  following  from  mesh 
connectivity.  Within  our  intrinsic  meshless  subdivision 
framework,  we  replace  these  extrinsic  weighted  aver¬ 
ages  of  topological  neighbours  by  intrinsic  weighted 
averages  of  intrinsically  neighbouring  points.  More 
specifically,  we  suggest  the  following  set  of  rules: 


Geometric  averaging  rule:  At  each  iteration,  replace 


Figure  4:  Updating  of  an  intrinsic  Voronoi  diagram  fol¬ 
lowing  insertion  of  a  point  pr  The  new  Voronoi  re¬ 
gion  R(pj,Pi)  is  grown  across  the  refined  point  cloud  P/ 
rendered  using  an  enlarged  point  size.  This  is  achieved 
by  propagating  a  front  from  the  newly  inserted  point  pj 
outwards,  (a)-(c),  until  it  encounters  its  neighbouring  re¬ 
gions  and  reaches  its  final  size,  (d). 


p  £  Pi  by  the  weighted  geodesic  centroid,  c(Np)  6%, 
of  its  intrinsic  neighbourhood  Np. 

Refinement  rule:  For  each  neighbour  qi  £  Np, 
consider  the  joint  intrinsic  neighbourhood  Npqi  of 
p.qi  £  P[.  Upsample  /)  by  inserting  the  weighted 
geodesic  centroid,  c(Npqj)  £  Pi+i,  of  Npq. . 


This  use  of  weighted  centroids  in  the  refinement  and 
geometric  averaging  rules  is  reminiscent  of  both  classi¬ 
cal  subdivision  schemes  [ZSOO]  (Section  2)  and  the  “re¬ 
peated  averaging”  approach  towards  the  generation  of 
subdivision  surfaces  ([OS03]  and  references  therein). 

By  performing  the  averaging  intrinsically  on  the  un¬ 
derlying  surface  represented  by  point  cloud  P,  the  above 
set  of  rules  raises  the  questions  of  how  to  compute 
geodesic  averages  on  manifold  surfaces  and  how  to  de¬ 
termine  a  suitable  weighting  scheme.  In  this  paper, 
we  are  interested  in  showing  the  conceptual  viability  of 
meshless  subdivision.  We  therefore  guide  the  choice  of 
weights  by  experimental  results  rather  than  any  theoret¬ 
ical  evidence  for  the  scheme’s  convergence  towards  a 
smooth  limit  surface.  Future  work  will  consider  formal 
proofs  of  the  scheme’s  convergence  to  a  limit  surface 
and,  consequently,  any  light  such  proofs  may  throw  on 
the  optimal  choice  of  weights. 

As  regards  the  computation  of  geodesic  centroids. 
Buss  and  Fillmore  [BF01]  present  an  algorithm  for  the 
computation  of  geodesic  averages  on  spheres.  We  gen¬ 
eralise  the  underlying,  earlier  idea  ([Kar77]  and  refer¬ 
ences  therein)  of  minimising  a  least  squares  expression 
in  weighted  geodesic  distances  to  non-spherical  geome¬ 
try  in  the  following  section. 
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Figure  5:  The  unweighted  centroid  of  a  (blue)  subset 
of  this  symmetric  set  of  points  on  a  surface  is  expected 
to  be  located  on  or  near  the  underlying  surface.  Since 
it  is  based  on  the  consideration  of  intrinsic  distances, 
this  is  the  case  when  computing  the  geodesic  centroid 
(red).  By  contrast,  in  the  case  of  the  Euclidean  averag¬ 
ing  of  the  (blue)  points,  the  resulting  centroid  (grey)  is 
located  away  from  the  underlying  surface.  This  effect 
gets  more  pronounced  when  increasing  the  size  of  the 
subset  whilst  the  position  of  the  geodesic  centroid  re¬ 
mains  unaffected  due  to  the  symmetry  of  the  point  con¬ 
stellation  (from  left  to  right). 


4.4  Geodesic  centroid  computation 

When  considering  the  concept  of  meshless  geometric 
subdivision  by  recursive  weighted  averaging  of  local 
point  neighbourhoods,  the  benefit  of  performing  these 
centroid  computations  intrinsically  rather  than  in  the 
Euclidean  domain  may  not  be  immediately  clear.  To 
see  the  need  for  operating  intrinsically  throughout,  note 
that  subdivision  should  ideally  generate  incresasingly 
smoother  representations  which  are  still  geometrically 
close  to  the  surface  represented  by  the  input  data.  This 
is  not  guaranteed  to  be  the  case  when  considering  Eu¬ 
clidean  instead  of  geodesic  centroids.  For  the  simple 
example  illustrated  in  Figure  5,  Euclidean  averaging 
moves  the  centroid  away  from  the  surface.  By  contrast, 
since,  by  definition,  geodesic  centroid  computation  is 
based  on  intrinsic  distances,  the  geodesic  centroid  falls 
onto  the  surface.  We  therefore  consider  geodesic  cen¬ 
troids  throughout  the  subdivision  process  and  present 
our  method  for  computing  these  centroids  on  manifold 
surfaces  next. 

The  weighted  (geodesic)  centroid  of  a  set  of  points 
{p i , . . .  ,p„}  CMC  M"'  is  defined  as  the  point  g  £  M 
which  minimizes 

1  " 

J(g)  ■=  xYwkdM(g,Pk), 

z  jfc=i 

where  w i, . . .  ,wn  are  the  weights:  0  <  w,-  <  1,  Ef=t  w,-  = 
1  and  is  the  geodesic  distance.  In  general,  the 

set  of  minimizers  of  /(•)  need  not  be  a  single  point. 
However,  when  pi,...,pn  are  all  contained  in  a  suffi¬ 
ciently  small  open  geodesic  ball  Bm  on  M,  there  exists 
exactly  one  minimizer,  ggM  of  J(-),  which  happens  to  lie 
in  Bm  [Kar77].  The  property  we  are  alluding  to  here  is 


convexity:  For  any  p.q  £  Bm,  we  require  that  the  short¬ 
est  geodesic  from  p  to  q  is  unique  in  M  and  contained  in 

BM- 

In  the  Euclidean  case,  direct  differentiation  of  J(-) 
yields  the  minimizer  g  =  ELi  wkPk-  This  simple  result 
does  not  extend  to  the  general  case  considered  here  but 
we  can  prove  that  any  minimizer  must  satisfy: 

V(g)  ■=  E  WkVM\dl[{g,pk)  =  0. 

k=  1  Z 


Then,  starting  from  a  good  initial  guess  go,  we  can  track 
the  minimizer  g  using  back  propagation  with  velocity 
field  F(-).  This  is  a  sensible  procedure  since  if  go  £  Bm 
and  Bm  as  above,  then  —V(x)  points  towards  ggM  for 
x  £  Bm  [Kar77]. 

In  practise,  we  set  go  =  II m  (ELi  WkPk),  where  YIm  '■ 
£1'm  — >  M  is  the  orthogonal  projection  operator.  We  now 
show  that  in  the  light  of  the  considerations  presented 
above,  this  represents  either  a  reasonable  initial  condi¬ 
tion  for  the  back  propagation  or  a  first  approximation  to 
the  intrinsic  centroid.  By  a  simple  application  of  Lemma 
17  of  [WD03],  we  have  that 


n 

(  n  \ 

Yj  WkPk  -  n M 

E  WkPk 

k=  1 

\k=\  ) 

<  C{diam(B))2 , 


where  C  is  a  global  constant  which  depends  on  the  cur¬ 
vatures  of  M.  Then,  let  B  =  Bm  (x,  £) ,  for  some  x  £  M 
and  e  >  0.  Since  || /?,-  —  x||  <  dM{pi,x)  <  £,  one  also 
has  \\Y!k=iwkPk~x\\  <  £•  Therefore,  since  ||go—  x||  < 
||g0  -  EL i  wkPk\ |  +  II  ELt  WkPk-x ||,  we  obtain 

||go  — x||  <Ce2  +  e, 


which  implies  dM{go,x)  <  e(1  +  _De)(1  +Ce),  for  an¬ 
other  constant  D  depending  on  global  metric  properties 
of  M  [MS03].  We  only  care  for  a  simplified  bound 


dM(go,x )  <£e. 

Finally,  let  8  >  0  be  the  maximal  8  >  0  such  that 
Bm(x, 8)  is  convex.  Note  that  it  is  a  fact  that  if  8  < 
^  min  (in j{M) ,  ^ J ,  where  in  j(M )  is  the  injectivity  ra¬ 
dius  of  M  and  K  bounds  all  sectional  curvatures  in  M 
from  above,  then  Bm(x.  8)  is  convex  for  any  x  £  M.  See 
§7.6  and  §7.7  in  [Cha97].  For  such  a  8  >  0  and  provided 
£  <  8 /E,  and  {p\,. . . ,pn }  C  Bm(x,e)  for  some  x  £  M, 
go  £  Bm(x,  8)  and  -V(go)  will  be  pointing  towards  gsM- 
Also,  in  case  we  want  to  use  go  as  an  approximation  to 
gsM,  we  have  the  (weak)  bound  dM(gBM,go)  <{E  +  1)e. 
Therefore,  go,  as  defined  above,  represents  a  sensible 
choice  as  the  initial  condition  of  an  eventual  back  prop¬ 
agation  step,  or,  in  any  case,  a  rough  approximation  to 
g bm  with  known  error  bound.  Note  in  particular  that  it 
is  also  a  useful  choice  from  the  point  of  view  of  compu¬ 
tational  ease. 
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To  demonstrate  the  applicability  of  this  approach  in 
the  context  of  meshless  subdivision,  we  consider  the 
case  of  M  representing  the  unit  sphere  in  the  following 
section. 


5  Experimental  results 

We  begin  by  considering  the  intrinsic  meshless  subdivi¬ 
sion  of  a  set  of  points  sampled  relatively  regularly  uni¬ 
formly  from  the  surface  of  the  unit  sphere.  This  ini¬ 
tial  restriction  to  spherical  geometry  allows  us  to  dis¬ 
cuss  qualitative  and  quantitative  aspects  of  our  meshless 
subdivision  operator  without  the  need  for  taking  into  ac¬ 
count  the  effects  of  any  particular  normal  estimation, 
projection  and  gradient  descent  operations  performed 
when  processing  more  complex  geometry. 

To  implement  the  intrinsic  centroid  computation 
method  and  thus  our  intrinsic  meshless  subdivision  op¬ 
erator  presented  in  the  previous  sections,  techniques  for 
the  computation  of  intrinsic  distances  between  points  on 
the  surface,  the  projection  of  the  starting  point  for  the 
back  propagation  onto  the  surface  and  the  computation 
of  the  back  propagation  itself  are  required.  In  the  case 
of  the  unit  sphere,  these  techniques  are  readily  avail¬ 
able.  Intrinsic  distances  between  points  follow  trigono¬ 
metrically  and  orthogonal  projection  is  trivial.  Simi¬ 
larly,  the  exponential  map  and  its  inverse  are  directly 
available  and  may  be  used  to  implement  the  back  prop¬ 
agation  procedure.  As  a  result,  for  the  case  of  spherical 
geometry,  our  approach  for  geodesic  centroid  computa¬ 
tion  narrows  down  to  the  technique  of  [BF01]. 

We  apply  our  intrinsic  meshless  subdivision  operator 
to  a  base  point  set  Pq  of  2144  points  sampled  relatively 
regularly  uniformly  from  the  unit  sphere  and  shown  in 
Figure  7.  To  obtain  initial  natural  neighbour  proximity 
information  for  points  in  Po,  we  use  the  intrinsic  point 
cloud  simplification  algorithm  [MD03]  without  request¬ 
ing  any  simplification.  The  application  of  our  intrinsic 
subdivision  operator  to  Pq  using  this  initial  proximity 
information  then  yields  the  subdivided  point  set  Pi  pre¬ 
sented  in  Figure  7.  The  result.  Pi,  obtained  from  the  ap¬ 
plication  of  the  operator  to  Pi  using  natural  neighbour 
information  updated  as  described  in  Section  4.2  is  also 
shown.  Given  the  relatively  strong  regularity  of  the  in¬ 
put  data,  uniform  weighting  was  used  for  both  the  re¬ 
finement  and  the  geometric  averaging  rule  in  both  itera¬ 
tions.  The  results  produced  by  our  meshless  subdivision 
operator  are  presented  alongside  the  point  sets  produced 
by  the  application  of  the  Loop  subdivision  scheme  to 
a  triangular  mesh  representation  of  the  base  point  set. 
As  indicated  by  the  detail  views  of  Figure  7,  the  point 
distributions  obtained  from  the  two  operators  after  two 
iterations  are  qualitatively  similar  with  the  slight  irregu¬ 
larities  in  the  distribution  of  Pq  being  slightly  more  pro¬ 
nounced  in  the  case  of  meshless  subdivision  due  to  the 
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Figure  6:  Histograms  of  the  (spherical)  distance  from 
each  point  in  Pi  (top)  and  Pi  (bottom)  to  its  closest 
neighbour  in  the  respective  set. 

use  of  uniform  weights.  The  size  of  8570  (Pi)  and  34275 
(P2)  points  respectively  of  the  subdivided  point  sets  co¬ 
incide  for  both  subdivision  schemes.  There  are  no  no¬ 
ticeable  differences  in  the  smoothing  effect  of  these  two 
operators. 

In  order  to  add  some  quantitative  detail  to  the  anal¬ 
ysis  of  the  point  sets  generated  by  our  subdivision  op¬ 
erator,  we  compute  the  mean  and  the  standard  devia¬ 
tion  of  the  distance  from  each  point  in  the  set  to  its 
closest  neighbor(s)  for  the  subdivided  point  sets  Pi  and 
P2.  For  each  x  in  the  point  set  P  C  S,  let  ,vt4  (x)  de¬ 
note  the  distance  from  x  to  its  kth  closest  neighbour. 
As  an  indicator  for  the  uniformity  of  the  density  of  a 
point  set  P  (at  level  k),  we  use,  for  a  positive  integer 
k,  p(k)  =  ■  Since  p(k)  represents  an  absolute 

measure  it  may  be  too  sensitive,  therefore  we  compute 
instead  p(k)  =  mean( '4)-std(Wt)  wjjere  mean  ancj  std 

stand  for  the  mean  and  standard  deviation  over  the  point 
set,  respectively. 

The  histograms  of  sd i(x)  corresponding  to  the  two 
sets  of  points  are  given  in  Figure  6.  Note  in  Table  1 
that  the  values  of  p(k)  (for  1  <  k  <  10)  are  quite  close 
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Model\  k 

1 

2 

3 

4 

5 

P\ 

0.807 

0.832 

0.832 

0.821 

0.788 

Pi 

0.781 

0.815 

0.823 

0.804 

0.785 

6 

7 

8 

9 

10 

Pi 

0.809 

0.8132 

0.786 

0.818 

0.828 

Pi 

0.799 

0.798 

0.784 

0.81 

0.822 

Table  1:  Values  of  p(k)  for  Pi  and  P2  and  with  k  £ 
{1,2,..., 10}. 

to  1  therefore  indicating  small  dispersion  up  to  the  10th 
closest  neighbor. 

Finally,  in  Figure  8  we  present  examples  dealing 
with  more  complicated  geometry.  These  results  were 
produced  by  using  the  projected  Euclidean  centroid  go 
(Section  4.4)  and  k  nearest  neighbourhoods.  It  is  easy 
to  prove  that  the  use  of  k  nearest  neighbourhoods  re¬ 
sults  in  non-uniform  point  distributions  even  if  starting 
from  a  uniform  sampling  such  as  the  one  produced  by 
the  technique  discussed  in  Section  4.1.  We  are  currently 
introducing  the  back  propagation  flow  to  produce  the 
true  geodesic  centroid  gsM  from  go,  and  are  changing 
the  neighbourhood  rule.  This  yields  a  more  uniform 
point  distribution  and  corresponding  examples  will  be 
reported  elsewhere.  Note  that  for  the  processing  of  the 
spherical  geometry,  we  use  ggM  and  the  neighbourhood 
is  completely  intrinsic. 

6  Conclusion 

In  this  paper,  we  introduced  the  concept  of  meshless, 
or  point  cloud,  subdivision.  Apart  from  the  introduc¬ 
tion  of  this  paradigm,  we  presented  a  particular  mesh¬ 
less  subdivision  scheme  based  on  the  computation  of 
intrinsic  centroids  for  point  clouds.  We  implemented 
and  showed  the  applicability  of  this  technique  with  the 
help  of  a  new  method  for  the  computation  of  geodesic 
centroids  on  manifold  surfaces. 

The  consideration  of  local  proximity  instead  of  mesh 
connectivity  information  makes  meshless  subdivision 
attractive  for  the  processing  of  point  clouds  represent¬ 
ing  higher-dimensional  surfaces.  By  operating  directly 
across  the  point  cloud,  problems  associated  with  the 
complexity  of  the  topological  subdivision  of  higher¬ 
dimensional  meshes  are  avoided.  Meshless  subdivision 
operators  may  be  devised  for  this  type  of  input  data 
which  upsample  the  point-sampled  geometry  more  con¬ 
servatively  than  the  operator  suggested  in  this  paper.  In 
this  direction,  we  are  considering  introducing  adaptive 
neighbourhoods  based  on  curvature  estimators  such  as 
those  reported  in  [CP03,  MN03],  We  leave  this  to  fu¬ 
ture  work. 

As  regards  future  work  on  the  theoretical  analysis 
of  our  meshless  intrinsic  subdivision  scheme,  Wallner 


and  Dyn  [WD03,  WD04]  show  the  convergence  and 
smoothness  of  non-linear  geodesic  curve  subdivision  by 
proximity  to  a  corresponding  linear  extrinsic  subdivi¬ 
sion  scheme.  Since  our  meshless  subdivision  scheme 
operates  inside  a  small  Euclidean  offset  band  around  the 
point  cloud  throughout,  following  [MS01,  MS03],  we 
intend  to  exploit  this  idea  in  a  surface  context  to  anal¬ 
yse  the  convergence  and  continuity  of  intrinsic  meshless 
subdivision  schemes. 
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Figure  7:  The  base  point  set,  Pq,  of  2144  points  acquired  from  the  unit  sphere  (left),  is  subdivided  recursively 
using  both  our  meshless  subdivision  operator  (bottom)  and  Loop  subdivision  (top).  The  triangular  base 
mesh  generated  from  Pq  for  the  support  of  Loop  subdivision  is  shown  on  the  left.  The  subdivided  point 
sets  P\  (8570  points)  and  Pi  (34275  points)  are  shown  in  the  center  and  on  the  right  respectively.  The 
corresponding  reconstructed  surfaces  are  given  next  to  each  point  set.  Both  these  surfaces  and  the  base 
mesh  were  computed  with  the  help  of  public  domain  software  [Par].  (This  figure  is  best  seen  on-screen.) 


Figure  8:  Example  for  the  meshless  geometric  subdivision  of  more  elaborate  geometry  in  form  of  the  Venus  data 
set.  The  base  point  set  is  displayed  on  the  left,  the  result  obtained  after  two  iterations  of  meshless 
geometric  subdivision  is  shown  on  the  right.  See  text  for  details  on  the  steps  for  the  production  of  these 
results  and  comments  on  the  slightly  non-uniform  distribution  of  the  subdivided  point  sets.  All  surfaces 
were  reconstructed  with  the  help  of  public  domain  software  [Par], 
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